The energetics of oxide surfaces by quantum Monte Carlo

Density functional theory (DFT) is widely used in surface science, but for some properties the predictions depend strongly on the approximation used for exchange–correlation energy. We note recent suggestions that the widely used generalized gradient approximation (GGA) is inferior to the local density approximation (LDA) for the surface formation energy σ of both transition metals and oxides. We report quantum Monte Carlo calculations of σ for the MgO(001) surface which support the accuracy of LDA for this case, and indicate that GGA is too low by ∼30%. We point out the potentially important implications of this result for nanoscience modelling. For many years, electronic-structure techniques have played a major role in surface science [1]. The most widely used of these techniques is density functional theory (DFT) [2], which is routinely employed to study the surfaces of metals, semiconductors and insulators. Recently, scattered pieces of evidence have appeared [3, 4] suggesting that commonly used DFT approximations may sometimes be seriously in error for basic quantities such as surface formation energies. Surprisingly, the evidence indicates that the local density approximation (LDA) of DFT, which is often an inferior approximation for other purposes, may in some cases be much better than the generalized gradient approximation (GGA) for the calculation of surface energies. Unfortunately, the paucity of reliable experimental data for surface energies makes it difficult to tell whether this is generally true. Quantum Monte Carlo [5] may offer a valuable way forward, since it is usually much more reliable than DFT, and can therefore be used to judge the accuracy of different DFT approximations. We report here QMC calculations of the MgO(001) surface formation energy, which support the superiority of LDA for this case. The surface formation energy σ of a material is the reversible work needed to cleave the bulk crystal, divided by the total area of new surface thus formed, the value of σ depending on 0953-8984/06/350435+06$30.00 © 2006 IOP Publishing Ltd Printed in the UK L435 L436 Letter to the Editor the crystallographic orientation of the surface. In DFT calculations, the standard procedure is to calculate the energy of a slab of material in periodic boundary conditions. With Eslab the total energy of the slab per repeating cell, and Ebulk the energy of the same number of atoms of bulk material, σ is then given by (Eslab − Ebulk)/A, where A is the total surface area (both faces) of the slab per repeating cell. DFT calculations of this kind on transition metals suggested many years ago [6] that LDA gives reasonably accurate σ values. A very recent re-examination of these calculations for six metals confirms this [8], and suggests that, roughly speaking, the widely used Perdew–Burke–Ernzerhof (PBE) form of the GGA [7] tends to underestimate σ by ∼30%. However, this conclusion can be challenged, because the ‘experimental’ values of σ are estimated indirectly from data on the surface tension of liquids. An important exception to this criticism is the case of Pb, where experimental σ values for three surfaces have been determined by scanning tunnelling microscopy (STM) measurements on Pb crystallites [9]. For this case, DFT calculations [3] show clearly that LDA reproduces the experimental σ data accurately, whereas GGA underestimates them by ∼30%. For oxides, it is also a familiar fact that GGA σ values are ∼30% lower than LDA [10], but the almost complete lack of experimental data makes it difficult to decide which is correct. One of the few exceptions is MgO, where the rather scattered data on the (001) surface [11] tend to favour LDA. Indirect confirmation for large GGA errors in oxide surface energies comes from a study [4] of the work of adhesion Wadh of Pd(111) to α-Al2O3(0001), for which accurate measurements are available. (Here, Wadh is the reversible work per unit area needed to separate the system containing the oxide– metal interface into its metal and oxide constituents.) The GGA value Wadh = 1.6 J m−2 is far below the LDA and experimental values of 2.4 and 2.8 J m−2. The authors argue [4] that the GGA error comes roughly equally from errors in the surface energies of the metal and the oxide. Fundamental insight into DFT errors for surface energetics comes from work on the jellium surface [12, 13]. Jellium is the homogeneous interacting electron gas neutralized by a uniform background; its density n is characterized by the mean inter-electron distance rs , defined by (4πr 3 s /3)n = 1, with rs in atomic units. The planar jellium surface is formed by having the neutralizing background occupy only the half-space x < 0, so that the electron number density n(x) in the ground state goes to its bulk value n for x → −∞ and to 0 for x → ∞. Accurate results for the formation energy of the jellium surface have been obtained [13] by extrapolating QMC calculations on neutral jellium spheres of different radius R. The extrapolation was performed by studying the large-R behaviour of the difference between the QMC total energy and the total energy calculated with DFT approximations, the main such approximations being the LDA [2], the GGA in the PBE form [7] and the meta-GGA [14]. This jellium surface work showed that: (i) for 2 rs 5, typical of simple metals such as Al and Na, the meta-GGA gives a very accurate surface energy σ , followed closely by the LDA, with the GGA being too low by a significant amount; (ii) as rs falls below 2, the GGA errors rapidly worsen. In the region rs ∼ 1.5, characteristic of transition metals and many oxides, the GGA σ is too low by ∼0.3 J m−2, a serious error because the surface energies of these materials are themselves in the region of 1 J m−2. It is convenient and instructive to present our QMC calculations on σ for MgO(001) by focusing on the difference between QMC and DFT values of the total energy of the system. The reason for doing this is that DFT errors for surface energetics are expected to be localized in the surface region [15]. This implies that thin slabs, containing only a few atomic layers, should suffice to assess the difference between the surface energy given by QMC and by DFT approximations, so that this difference will converge rapidly with increasing slab thickness. In general, the surface energy of a material is significantly influenced by relaxation of ions in the surface region. It is straightforward to include this relaxation in DFT calculations, and this has Letter to the Editor L437 been done in the present work. However, structural relaxation cannot be routinely performed in the QMC framework, and for our QMC calculations we have taken the relaxed structure of the slab from our DFT calculations. Fortunately, in the case of MgO(001), relaxation is very small, and has a negligible effect on σ [16], so this procedure should be completely satisfactory. In the following, we report first a set of DFT calculations of σ , stressing the technical issues that must be brought under control to achieve reliable results. Then we describe the QMC calculations themselves, noting what must be done to ensure reliable results for the difference of DFT and QMC σ values. Our DFT calculations used the standard pseudopotential/plane-wave technique [2], and were performed using the VASP code [17]. The surface was modelled using periodically repeated slab geometry, the calculation conditions being characterized by basis-set completeness (plane-wave cut-off energy Ecut), Brillouin-zone (BZ) sampling of the electronic states, the width L of the vacuum layer separating successive slabs and the number of ionic layers Nlayer in each slab. The formula σ = (Eslab − Ebulk)/A applies for all Nlayer 1. The bulk energy Ebulk is Nlayer times the bulk energy per layer ebulk, and it is convenient to obtain ebulk from the difference of Eslab values for successive values of Nlayer in the limit of large Nlayer. For given Nlayer, we always insist on convergence of the calculated σ with respect to Ecut, BZ sampling and L to within 0.01 J m−2 (this tolerance is satisfied for L > 6 Å). As expected from earlier work, σ for MgO(001) converges rapidly with respect to Nlayer, the residual errors being below 0.01 J m−2 for Nlayer 2. For a given DFT approximation, the calculated σ depends a little on MgO lattice parameter a0. For the experimental value a0 = 4.21 Å, we obtain σ = 1.24 and 0.87 J m−2 with LDA and GGA(PBE) respectively. The finding that the GGA value is ∼30% lower than the LDA one is typical of oxide surfaces. The calculation of σ by QMC is not standard, and we are not aware of previous calculations of σ for any oxide surface using QMC, though our recent QMC calculations on perfect and defective MgO crystals [18] indicated the feasibility of the present calculations. We refer the reader to reviews for details of QMC (e.g. [5]). We recall that for high-precision results it is essential to use diffusion Monte Carlo (DMC), in which the many-electron wavefunction is evolved in imaginary time, starting from an optimized trial wavefunction generated in prior variational Monte Carlo calculations. The only error inherent in DMC is ‘fixed-node’ error, due to the fact that the nodes of the many-electron wavefunction are constrained to be those of the trial wavefunction. For many systems, including jellium, the evidence is that the fixednode error is extremely small. For wide-gap systems such as MgO, the errors should be no greater. Our calculations were performed with the CASINO code [19], using the same Hartree–Fock pseudopotentials as in our previous work [18]. The trial wavefunctions were of the usual Slater–Jastrow type, with single-electron orbitals obtained with the plane-wave code PWSCF [20], generally using the large plane-wave cut-off of 4082 eV. These orbitals were represented in CASINO using the recently reported ‘blip-function’ real-space basis set [21]. The DMC calculations all used a time-step of 0.005 au, and a mean number of walkers equal to 10 240. The calculations were done with free boundary conditions (i.e. no periodicity) normal to the surface. Our DMC calculations were performed on a series of MgO slabs with the number of layers Nlayer running from 1 to 5. For each Nlayer, convergence must be demonstrated with respect to basis-set completeness and the size of the repeating surface unit cell. Basis-set errors with the blip basis set are readily made negligible, as shown earlier [21]. In DMC calculations, the wavefunctions are real, so that BZ sampling is generally impossible, and calculations are usually performed at the -point; this is why convergence with respect to size of the surface unit cell must be checked. Our main DMC calculations used the 2 × 2 surface unit cell, for L438 Letter to the Editor which the repeating cells contain from 16 (Nlayer = 1) to 80 (Nlayer = 5) ions; we show below that larger surface cells would give almost identical results. The raw output from these calculations is DMC total energies EDMC slab for the five Nlayer values. Following our strategy, we now study the difference Eslab ≡ E slab − EDFT slab , with the DFT slab energy calculated with exactly the same slab and the same ( -point) BZ sampling as in the DMC calculations. Since the jellium results indicate that LDA surface energies are likely to be closer to DMC than those from GGA, we use LDA values for EDFT slab . When plotted against Nlayer, Eslab will tend asymptotically (Nlayer → ∞) to a straight line, whose slope is equal to the difference between the DMC and LDA bulk energies ebulk per layer, and whose Nlayer = 0 intercept divided by A gives the difference of DMC and DFT surface energies σ ≡ σDMC − σDFT. Since ebulk is large, and since Eslab contains the statistical errors of DMC, it is helpful to start this analysis by performing a least-squares straight-line fit a + bNlayer to the values of Eslab, and then to use the resulting b value to form the quantity ̃Eslab ≡ Eslab − bNlayer. The Nlayer → ∞ straight-line asymptote of ̃Eslab has the same Nlayer = 0 intercept as that of Eslab. For Nlayer = 1, 2, . . . , 5, we find the five values ̃Eslab = −0.019(2), −0.009(6), −0.007(9), −0.011(13) and −0.014(15) J m−2. This immediately shows that the DMC and LDA values of σ are almost exactly the same. Within our rather small statistical errors of at worst 0.015 J m−2, the difference between the DMC and LDA surface energies has the very small value of −0.01 J m−2. To check the errors due to use of the 2 × 2 surface cell (i.e. errors of BZ sampling), we have performed LDA calculations on slabs having large surface cells with a series of Nlayer values, using -point sampling. The σ values thus obtained differ from the LDA σ value extracted by similar -point calculations with the 2×2 surface cell by only 0.01 J m−2. The LDA value of σ for the lattice parameter we are using, converged with respect to BZ sampling and slab thickness, is 1.20 J m−2, and we conclude from our ̃Eslab values that the fully converged DMC value is 1.19± 0.01 J m−2. MgO is one of the few oxides for which experimental values of the surface energy are available [11]. Exploiting the fact that MgO cleaves readily along the (001) plane, the experiments measure the work of cleavage, thus ensuring that the results for σ cannot be influenced by surface contamination. Our σDMC value of 1.19 J m−2 is consistent with the measured values [11], which fall in the range 1.04–1.20 J m−2. Our QMC calculations thus support the accuracy of LDA for the surface formation energy σ of MgO(001), and indicate that GGA underestimates σ by ∼30%. This adds to the existing evidence for an approximate rule of thumb that LDA predictions of σ tend to be more accurate than GGA for quite a wide range of materials, including in particular transition metals and oxides. The present QMC calculations could rather easily be applied to a number of other oxides, including alumina and silica, and we hope to do this. In considering this possibility, it will be noted that it remains difficult to calculate forces within QMC, though there are signs that this problem may be soluble [22]. This means that it is difficult to calculate the relaxed structures of oxide surfaces directly from QMC. However, although σ for oxides is sometimes strongly influenced by relaxation, it is possible that the difference of DFT and QMC values may be less strongly influenced. It is interesting also to note another possibility for the future. Although we have focused here on QMC calculations, it is likely that post-Hartree– Fock calculations (specifically, at the MP2 level) will soon be feasible in periodic boundary conditions [23]. This will offer a completely independent way of testing DFT predictions of σ for oxides. We point out some important implications of our results. The energetics of nanocrystalline systems is strongly affected by surface energetics. One example is that nanocrystals can adopt crystal structures that are not the most stable bulk structures, because they are stabilized Letter to the Editor L439 by having low surface energies [24]. Now GGA is generally better than LDA for bulk energetics. If, therefore, the reverse is true for surface energetics, then DFT modelling of nanocrystals could face serious problems. The evidence suggests that this might be the case. Another example is the formation of nanostructures on surfaces. Sometimes, the energetics of nanostructures involves a competition between strain energy and surface energy [25]. If the accuracies of LDA and GGA are similar for strain energy but substantially different for surface energies then their predictions for the stability of nanostructures could differ greatly. We suggest that this cannot be ruled out. In summary, we have presented a set of QMC calculations of the surface formation energy σ of MgO(001), which are accurately converged with respect to the basis set and other technical parameters. The results have been used to test DFT predictions, and we confirm existing evidence that LDA values of σ are more accurate than those of GGA, which underestimates σ typically by ∼30%. We have noted some important implications for the DFT modelling of the energetics of nanostructures. The work was conducted as part of a EURYI scheme award as provided by EPSRC (see www. esf.org/euryi). Allocation of computer time at the HPCx national service was provided by the Materials Chemistry Consortium, the UKCP Consortium and the Mineral Physics Consortium. We thank the Texas Advanced Computing Center at the University of Texas at Austin for providing high-performance computing resources. MJG acknowledges useful discussions with M Scheffler.